/*  -*- mode: c++; coding: utf-8; c-file-style: "stroustrup"; -*-

    Copyright 2008 Asier Aguirre <asier.aguirre@gmail.com>
    This file is part of Intromit.

    Intromit is free software: you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation, either version 3 of the License, or
    (at your option) any later version.

    Intromit is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with Intromit.  If not, see <http://www.gnu.org/licenses/>.
*/

#ifndef MATRIX_HH
#define MATRIX_HH
//! \addtogroup introsight_framework
//@{

#include "util.hh"
#include "refcnt.hh"

class matrix: public reference_counter
{
    // temporal memory for minor matrices
    double *minor;
    int n_minor, max_minor;

    // STATUS VARIABLES
    int m/*rows*/, n/*columns*/, mxn;
    // matrix coeficients cij=c[i+j*m]
    double *c;

    // PRIVATE FUNCTIONS
    void ctor();
    void dtor();

    // get more memory for minor matrixs
    void minimum_minor_size(int s);
    // uses m [sxs] coefs and sets minor in minor[n_minor...]
    // it supposes minor has enough space at its end (do minimum_minor_size externally)
    // returns next minor free position
    int minor_matrix(double *m, int s, int i, int j);
    double cofactor_recur(double *m, int s, int ii, int jj);

public:
    matrix();
    virtual ~matrix();
    int init(int m/*rows*/, int n/*columns*/);

    // set to identity matrix
    void identity();

    // set/get matrix from its coefs (should be at least n*m)
    void set_coefs(double *coefs);
    inline void set_coef(int i, int j, double coef) {c[i+j*m]=coef;}
    inline double *get_coefs() {return c;} // volatile pointer
    void get_coefs_copy(double *c);
    inline double get_coef(int i, int j) {return c[i+j*m];}

    // vector operations
    // factor of size n; and res of size m
    void vect_mult(double *factor, double *res);
    // for square matrixs: res = this * m
    bool matrix_mult(matrix *m, matrix *res);
    bool matrix_mult(double *coefs, double *res); // coefs and res should be different

    // determinant
    double determinant();

    // cofactor
    double cofactor(int i, int j);

    // matrix inversion
    bool invert(matrix *result);
    bool invert(double *result);

    // print matrix to file
    void dump(FILE *arch=stdout, int colours=1, char *tab=(char*)"");
    inline void dump_stdout() {dump();}
    void dump_octave_stdout();
};

class matrix4x4: public reference_counter
{
    // STATUS VARIABLES
    // matrix coeficients cij=c[i+(j<<2)]; where i->row and j->column
    double c[16];

    // PRIVATE FUNCTIONS
    void ctor();
    void dtor();

    // transformation matrixs coefs
    static void rotation_coefs(double angle, double vector[3], double coefs[16]); // angle in rads
    static void scale_coefs(double xyz, double coefs[16]);
    static void scale_coefs(double xyz[3], double coefs[16]);
    static void translate_coefs(double orig[3], double coefs[16]);

    // uses m [sxs] coefs and sets minor in minor[n_minor...]
    // it supposes minor has enough space at its end (do minimum_minor_size externally)
    // returns next minor free position
    static void minor_matrix(double *m, double *result, int s, int i, int j);
    static double cofactor_recur(double *m, int s, int ii, int jj);

public:
    matrix4x4();
    matrix4x4(double coefs[16]);
    virtual ~matrix4x4();

    // set to identity matrix
    void identity();

    // set/get matrix from its coefs (should be at least n*m)
    void set_coefs(double coefs[16]);
    inline void set_coef(int i, int j, double coef) {c[i+(j<<2)]=coef;}
    inline double *get_coefs() {return c;}
    void get_coefs_copy(double c[16]);
    inline double get_coef(int i, int j) {return c[i+(j<<2)];}

    // vector operations
    // factor of size n; and res of size m (m=n=4)
    void vect_mult(double factor[4], double res[4]);
    static void vect_mult(double matrix_factor[16], double vect_factor[4], double res[4]);
    static void vect_mult_3x3_1(double matrix_factor[16], double vect_factor[3], double res[3]);
    static void vect_mult_3x3_0(double matrix_factor[16], double vect_factor[3], double res[3]);
    static void vect_mult_3x3_1_inv(double vect_factor[3], double matrix_factor[16], double res[3]);
    static void vect_mult_3x3_0_inv(double vect_factor[3], double matrix_factor[16], double res[3]);
    // for square matrixs: res = this * m
    void matrix_mult(matrix4x4 *m, matrix4x4 *res);
    void matrix_mult(double coefs[16], double res[16]); // coefs and res should be different
    void matrix_mult_inv_order(double coefs[16], double res[16]); // insted of this * coefs -> coefs * this

    // determinant
    double determinant();
    static double determinant(double coefs[16]);

    // cofactor
    double cofactor(int i, int j);
    static double cofactor(int i, int j, double coefs[16]);

    // matrix inversion
    bool invert(matrix4x4 *result);
    bool invert(double result[16]);
    static bool invert(double coefs[16], double result[16]);

    // matrix spatial typical matrixs
    void rotation(double angle, double vector[3]); // angle in rads
    void scale(double xyz);
    void scale(double xyz[3]);
    void translate(double orig[3]);

    // matrix spatial typical operations
    void do_rotation(double angle, double vector[3]); // angle in rads
    void do_rotation(double angle, double x, double y, double z); // angle in rads
    void do_scale(double xyz);
    void do_scale(double xyz[3]);
    void do_translate(double orig[3]);
    void do_translate(double x, double y, double z);

    // print matrix to file
    void dump(FILE *arch=stdout, int colours=1, char *tab=(char*)"");
    static void dump_coefs(double c[16]);
    void dump_stdout();
    void dump_octave_stdout();
};

//@}
#endif

/* Offline metadata attributes:

   @metadata_attrs(
	matrix
	{
		# members
		n, m, c:
						select, read_only

		# methods
		:
						select
	}

	matrix4x4
	{
		# members
		c:
						select

		# methods
		identity, set_coefs, get_coef,
		matrix_mult,
		matrix_mult_inv_order,
		determinant, cofactor, invert,
		rotation, scale, translate,
		do_rotation, do_scale,
		do_translate, dump_coefs,
		dump_stdout,
		dump_octave_stdout:
						select
	}
   )
*/

